#include "ParticleTracking.h"

ParticleTracking::ParticleTracking()
{
	tAll = make_shared<vector<double>>(1, 0.0);
	pVAll = make_shared<vector<double>>(1, 0.0);
	pXAll = make_shared<vector<double>>(1, 0.0);
}

ParticleTracking::ParticleTracking(InitialConditions pConditions, string outputFileName)
{
	tAll = make_shared<vector<double>>(pConditions.nSteps + 1, 0.0);
	pVAll = make_shared<vector<double>>(pConditions.nSteps + 1, 0.0);
	pXAll = make_shared<vector<double>>(pConditions.nSteps + 1, 0.0);

	(*pVAll)[0] = pConditions.pV0; (*pXAll)[0] = pConditions.pX0;

	for (int i = 0; i < pConditions.nSteps; i++)
	{
		int i_store = i + 1;
		(*tAll)[i_store] = (*tAll)[i] + pConditions.dt;
		double pReTemp = CalculateRe(pConditions.rhoG, pConditions.pDia, abs(pConditions.velG - (*pVAll)[i]), pConditions.muG);
		double CdTemp = CalculateCd(pReTemp);
		double pTaoTemp = CalculateTao(pConditions.pRho, pConditions.pDia, pConditions.muG, CdTemp, pReTemp);
		//Update particle velocity
		(*pVAll)[i_store] = (*pVAll)[i] + pConditions.dt * (pConditions.velG - (*pVAll)[i]) / pTaoTemp;
		(*pXAll)[i_store] = (*pXAll)[i] + pConditions.dt * (*pVAll)[i];
	}

	string outputFileDirect = GetOutputDirectory(outputFileName);
	ofstream outputFile(outputFileDirect);
	assert(outputFile.is_open());
	for (int i = 0; i < pConditions.nSteps + 1; i++)
	{
		if (i == 0)
		{
			outputFile << "time (s) " << "Vp (m/s) " << "X (m) " << endl;
		}
		outputFile << (*tAll)[i] << " " << (*pVAll)[i] << " " << (*pXAll)[i] << endl;
	}
	outputFile.close();
}

string ParticleTracking::GetOutputDirectory(string fileName)
{
	char homeDirect[1024];
	_getcwd(homeDirect, 1024);
	string localDirect = (string)(homeDirect)+(string)("\\LocalResults\\") + fileName;
	return localDirect;
}
double ParticleTracking::CalculateRe(double rhoG, double DiaP, double slipVel, double muG)
{
	return rhoG * DiaP * slipVel / muG;
}
double ParticleTracking::CalculateCd(double Re)
{
	// empirical spherical model by Mori and Alexander
	double Cd = 0.0;
	if (Re < 0.1)
	{
		Cd = 24.0 / Re;
	}
	else if (Re >= 0.1 && Re < 1.0)
	{
		Cd = 22.73 / Re + 0.0903 / (Re * Re) + 3.69;
	}
	else if (Re >= 1.0 && Re < 10.0)
	{
		Cd = 29.1667 / Re - 3.8889 / (Re * Re) + 1.222;
	}
	else if (Re >= 10.0 && Re < 100.0)
	{
		Cd = 46.5 / Re - 116.67 / (Re * Re) + 0.6167;
	}
	else if (Re >= 100.0 && Re < 1000.0)
	{
		Cd = 98.33 / Re - 2778.0 / (Re * Re) + 0.3644;
	}
	else if (Re >= 1000.0 && Re < 5000.0)
	{
		Cd = 148.62 / Re - (4.75e4) / (Re * Re) + 0.357;
	}
	else if (Re >= 5000.0 && Re < 10000.0)
	{
		Cd = -490.546 / Re + (57.87e4) / (Re * Re) + 0.46;
	}
	else if (Re >= 10000.0 && Re < 50000.0)
	{
		Cd = -1662.5 / Re + (5.4167e6) / (Re * Re) + 0.5191;
	}
	else
	{
		assert(false && "Cd value is out of Range!");
	}
	return Cd;
}
double ParticleTracking::CalculateTao(double rhoP, double DiaP, double muG, double Cd, double Re)
{
	return rhoP * pow(DiaP, 2) * 24.0 / (18.0 * muG * Cd * Re);
}
